Bienvenid@s a la segunda tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en inferencia estadística, diseño de experimentos y test de hipótesis. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Slides de las clases:
Enlaces a videos de las clases:
Documentación:
Usted deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:
Esta pregunta tiene como objetivo comprobar a través de gráficos las características que poseen los estimadores.Por favor responda de forma separada las siguientes preguntas:
Respuesta: En primera instancia, realizamos 50.000 experimentos aleatorios para una distribución de Bernoulli. Luego, graficamos como se va viendo el estimador \(\hat{p}_{n}\) mientras se va avanzando en la cantidad de experimentos.
n_max <- 50000
p <- 0.5
# Vectores del estimador y de la distribución de Bernoulli
estimator <- vector(length = n_max)
bernoulli <- rbinom(n_max, size = 1, prob = p)
# Para cada valor del vector, se promedia con los anteriores
# y se agregan al estimador
for (k in 1:n_max) {
estimator[k] <- sum(bernoulli[1:k]) / k
}
# Gráfico de los experimentos
plot_est <- plot_ly(
x = 1:n_max,
y = estimator,
mode = "lines",
type = "scatter",
name = "Experimentos",
line = list(color = "blue")
) %>%
add_trace(
x = 1:n_max,
y = 0.5,
name = "p = 0.5",
line = list(color = "red")
) %>%
layout(
title = "Estimador Pn para X ~ Bernoulli(p) con p = 0.5",
xaxis = list(title = sprintf("Pasos (%d)", n_max)),
yaxis = list(title = "Probabilidad")
)
plot_est
Podemos notar que a medida que vamos avanzando en la cantidad de experimentos, el estimador se acerca cada vez más al valor teórico de \(p = 0.5\), lo cual verifica que el estimador es consistente.
Por otra parte, realizamos el mismo procedimiento para el estimador \(T_{n}\).
n_max <- 50000
estimator <- vector(length = n_max)
bernoulli <- rbinom(n_max, size = 1, prob = 0.5)
noise <- rnorm(n_max)
T <- vector(length = n_max)
for (k in 1:n_max) {
estimator[k] <- sum(bernoulli[1:k]) / k
T[k] = estimator[k] + noise[k]
}
plot_est <- plot_ly(
x = 1:n_max,
y = T,
mode = "lines",
type = "scatter",
name = "Experimentos",
line = list(color = "blue")
) %>%
add_trace(
x = 1:n_max,
y = 0.5,
name = "p = 0.5",
line = list(color = "red")
) %>%
layout(
title = "Estimador Tn para X ~ Bernoulli(p) con p = 0.5",
xaxis = list(title = sprintf("Pasos (%d)", n_max)),
yaxis = list(title = "Probabilidad")
)
plot_est
Se puede apreciar que el estimador \(T_n\) no es consistente ya que no converge al valor esperado de \(p = 0.5\). La inconsistencia de este estimador se debe a que al estimador \(\hat{p}_n\) se le agrega una variable normal en cada iteración, un “ruido”.
plot_est <- plot_ly(
x = 1:n_max,
y = noise,
mode = "lines",
type = "scatter",
name = "Ruido",
line = list(color = "blue")
) %>%
layout(
title = "Ruido e ~ N(0, 1)",
xaxis = list(title = sprintf("Pasos (%d)", n_max)),
yaxis = list(title = "Probabilidad")
)
plot_est
Este ruido es el que no permite que el estimador \(T_n\) converga al parámetro \(p\), ya que hace que el estimador “oscile” aleatoriamente, evitando que converga y por lo tanto, haciendo de \(T_n\) un estimador insesgado pero inconsistente.
El objetivo de esta pregunta es visualizar los intervalos de
confianza en datos simulados de una población, para visualizar la
incertidumbre que presenta una estimación. Para esto, ustedes deberán
generar datos de una distribución exponencial, la cual deberán
considerar como los datos de la población. En base a los datos
generados, estimen la distribución de la media de la población a través
de la sampling distribution de la media. Notar que el valor
obtenido en cada muestra les entregará un estimador de la media, o sea,
para cada valor podremos calcular un intervalo de confianza. Hecho esto,
calculen el intervalo de confianza del \(95\%\) para cada una de las medias
estimadas, utilizando la función de cuartil vista en clases.
Para la elaboración de esta parte de la tarea, se recomienda realizar el experimento con la siguiente secuencialidad:
Notar: Responder cada una de las preguntas señaladas en esta sección.
Hints:
sampling distribution podría serle
útil el comando sample().
Del gráfico es posible observar que la línea punteada es la media de la población y los puntos de colores son las estimaciones con sus respectivos intervalos de confianza. Notar que para el plot no se utilizaron las 5000 veces, se recomienda utilizar 100 valores para visualizar bien el fenómeno.
Respuesta:
# Definimos tamaños de muestreo
tamano_muestra = 30
n_muestras = 5000
# Generamos una exponencial para luego generar el subsampling de ella
exponencial = rexp(10000, rate = 2)
# Obtenemos la media poblacional de la exponencial
# Sampling distribution, calculo del intervalo de confianza y proporción.
# - En base a la distribución exponencial, generamos multiples sampling
# distribution.
# - Se estima la media del muestreos y obtenemos el intervalo de confianza de
# cada una de las muestras
# Plot de Intervalos de confianza (ver respuesta esperada)
# Plot de proporción de Intervalos de confianza
En esta pregunta deberán trabajar con el dataset
Body Measurements_original.csv. El objetivo será visualizar
e inferir los parámetros que componen a dos variables del dataset. Para
esto deberá visualizar el comportamiento de la likelihood, utilizando
diferentes cantidades de datos, y realizar la optimización de la
likelihood para obtener los estimadores de las variables a través de la
función nlminb(). Notar que esta pregunta consiste en dos
partes.
TotalHeight, donde deberá asumir y realizar los
siguientes puntos:Age, y estimar a
traves de la -log(likelihood) solo los parámetros de la
distribución que observa (notar que solo debe inferir los estimadores de
variable escogida). Para señalar la distribución de los datos se
recomienda realizar el plot de densidad y comparar con el comportamiento
de las distribuciones teóricas vistas en clases.Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo, la librería MASS).
Respuesta
Primero, visualicemos la densidad de la variable
TotalHeight
data_raw <- read.csv("Body Measurements _ original.csv")
totalheight_plot <- ggplot(
data = data_raw,
aes(x = TotalHeight, color = "density")
)
totalheight_dens <- geom_density(color = "#67B7D1")
totalheight_hist <- geom_histogram(
aes(y = after_stat(density)),
bins = 10,
fill = "#67B7D1",
alpha = 0.5
)
final_plot <- totalheight_plot + totalheight_dens + totalheight_hist +
ylab("") + xlab("") + scale_color_manual(values = c("density" = "#67B7D1"))
ggplotly(final_plot) %>%
layout(
xaxis = list(title = "TotalHeight"),
yaxis = list(title = "Frequency")
)
Claramente se aprecia que esta variable se asemeja mucho a una distribución normal. Ahora, veamos un gráfico de calor para visualizar el comportamiento de \(\mu\) y \(\sigma\).
likelihood <- function(mu, sigma, sub = NULL) {
x_org <- data_raw$TotalHeight
if (is.null(sub))
x <- x_org
else
x <- sample(x_org, sub)
n <- length(x)
x_mean <- mean(x)
x_var <- var(x)
first_sum <- (n * x_var) / (2 * sigma^2)
second_sum <- (n * (x_mean - mu)^2) / (2 * sigma^2)
-n * log(sigma) - first_sum - second_sum
}
mu <- seq(from = 20, to = 80, by = 0.5)
sigma <- seq(from = 5, to = 23, by = 0.5)
ll_plot <- Vectorize(likelihood)
ll_plot <- outer(X = mu, Y = sigma, ll_plot)
filled.contour(x = mu, y = sigma, z = ll_plot,
xlab = expression(mu), ylab = expression(sigma))
Se puede apreciar que, muy probablemente, se tiene \(\mu \in (40, 60)\) y \(\sigma \in (5, 10)\).
Ahora, fijemos \(\sigma = 12\) y variemos el valor de \(\mu\)
sigma <- 12
p <- plot_ly(x = mu, y = likelihood(mu, sigma, 100),
type = "scatter", name = "100 datos")
p <- add_trace(p, y = likelihood(mu, sigma, 300), name = "300 datos")
p <- add_trace(p, y = likelihood(mu, sigma), name = "Todos los datos")
p <- p %>% layout(title = "Likelihood para valores de mu y sigma = 12",
xaxis = list(title = "mu"),
yaxis = list(title = "valor de likelihood")
)
p
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Se puede apreciar que el máximo valor de la función likelihood disminuye mientras aumenta la cantidad de datos considerada. No obstante, esto no influye en el óptimo para \(\mu\), ya que este sigue siendo el mismo para cualquier cantidad de datos.
Minimicemos la función likelihood para encontrar los valores de \(\mu\) y \(\sigma\).
# Primero, redefinimos la función likelihood para aplicar nlminb
likelihood <- function(param) {
mu <- param[1]
sigma <- param[2]
x <- data_raw$TotalHeight
n <- length(x)
x_mean <- mean(x)
x_var <- var(x)
first_sum <- (n * x_var) / (2 * sigma^2)
second_sum <- (n * (x_mean - mu)^2) / (2 * sigma^2)
n * log(sigma) + first_sum + second_sum
}
# Aplicamos nlminb
MLE <- nlminb(
objective = likelihood,
start = c(10, 10),
lower = c(0.001, 0.001),
upper = c(500, 500)
)
cat("mu = ", MLE$par[1], ", sigma = ", MLE$par[2])
## mu = 48.1187 , sigma = 12.15672
De aquí obtenemos que los valores son \(\mu = 48.118\) y \(\sigma = 12.156\).
Ahora, trabajaremos con la variable ShoulderToWaist.
Veamos su densidad:
data_plot <- ggplot(
data = data_raw,
aes(x = ShoulderToWaist, color = "density")
)
data_dens <- geom_density(color = "#67B7D1")
data_hist <- geom_histogram(
aes(y = after_stat(density)),
bins = 10,
fill = "#67B7D1",
alpha = 0.5
)
final_plot <- data_plot + data_dens + data_hist +
ylab("") + xlab("") + scale_color_manual(values = c("density" = "#67B7D1"))
ggplotly(final_plot) %>%
layout(
xaxis = list(title = "ShoulderToWaist"),
yaxis = list(title = "Frequency")
)
Notemos que esta variable también se asemeja a una distribución normal. Usando el mismo método, se ve el mapa de calor y se calculan los valores óptimos.
likelihood <- function(mu, sigma) {
x <- data_raw$ShoulderToWaist
sum(log(dnorm(x, mean = mu, sd = sigma)))
}
mu <- seq(from = 10, to = 30, by = 0.5)
sigma <- seq(from = 3, to = 10, by = 0.5)
ll_plot <- Vectorize(likelihood)
ll_plot <- outer(X = mu, Y = sigma, ll_plot)
filled.contour(x = mu, y = sigma, z = ll_plot,
xlab = expression(mu), ylab = expression(sigma))
# Se redefine la función likelihood para usar nlminb
likelihood <- function(param) {
mu <- param[1]
sigma <- param[2]
x <- data_raw$ShoulderToWaist
-sum(log(dnorm(x, mean = mu, sd = sigma)))
}
MLE <- nlminb(objective = likelihood, start = c(10, 10), lower = c(0.001, 0.001), upper = c(500, 500))
cat("mu = ", MLE$par[1], ", sigma = ", MLE$par[2])
## mu = 17.90084 , sigma = 5.375555
De aquí obtenemos que los valores encontrados son \(\mu = 17.9\) y \(\sigma = 5.375\).
Para esta pregunta será necesario cargar el dataset
SAT_GPA.csv, con el que estudiaremos la correlación entre
las variables SAT y GPA. Dentro de las variables: GPA representa el
rendimiento académico de un estudiante en el sistema estadounidense,
mientras que SAT es una prueba de admisión universitaria en estados
unidos. Las actividades por realizar en esta pregunta son las
siguientes:
Nota: No se permite la utilización de librerías de bootstrap para esta parte.
Respuesta:
A work by CC6104